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Abstract 

This paper presents the results of a study on the thermomechanical behavior of internally cooled 
silicon nitride structures. Silicon nitride is under consideration for elevated temperature aerospace 
engine applications, and techniques for lowering the operating temperature of structures 
composed of this material are under development. Lowering the operating temperature provides a 
large payoff in terms of fatigue life and may be accomplished through the use of thermal barrier 
coatings (TBC’s) and the novel concept of included cooling channels. Herein, an in-depth study is 
performed on the behavior of a flame-impinged silicon nitride plate with a TBC and internal 
channels cooled by forced air. The analysis is performed using the higher order theory for 
functionally graded materials (HOTFGM), which has been developed through NASA Glenn 
Research Center funding over the past several years. HOTFGM was chosen over the traditional 
finite element approach as a prelude to an examination of functionally graded silicon nitride 
structures for which HOTFGM is ideally suited. To accommodate the analysis requirements of 
the internally cooled plate problem, two crucial enhancements were made to the two-dimensional 
Cartesian-based version of HOTFGM, namely, incorporation of internal boundary capabilities 
and incorporation of convective boundary conditions. Results indicate the viability and large 
benefits of cooling the plate via forced air through cooling channels. Furthermore, cooling can 
positively impact the stress and displacement fields present in the plate, yielding an additional 
payoff in terms of fatigue life. Finally, a spin-off capability resulted from inclusion of internal 
boundaries within HOTFGM; the ability to simulate the thermo-elastic response of structures 
with curved surfaces. This new capability is demonstrated, and through comparison with an 
analytical solution, shown to be viable and accurate. 


1 .0 Introduction 

Cooled ceramics are under consideration for high temperature aircraft engine applications. 
Ceramics offer higher operating temperatures than metals, which are the traditional choices for these 
applications, but due to ceramics’ low thermal conductivity, thermal shock induced cracking can be a 
problem. By lowering the operating temperature of ceramic components, through both thermal barrier 
coatings (TBC’s) and internal cooling, the driving force for thermal shock can be reduced while the 
resistance to cracking and fatigue failure can be increased dramatically. One particular ceramic that shows 
promise as a turbine blade material is silicon nitride (SLN 4 ). Figure 1 provides sample-applied stress vs. 
life curves for SLN 4 . As Fig. 1 shows, a large payoff in fatigue life can be realized by lowering the 
operating temperature, even by as little as two hundred degrees F [1], 
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The concept of employing a low thermal conductivity TBC on a hot surface to lower the 
operating temperature of engine components is well known. Recent work has indicated that including 
internal cooling channels within ceramic engine components and forcing air through the channels might 
further reduce component-operating temperatures [2]. A study to evaluate the feasibility and efficacy of 
internally cooled Si 3 N 4 is in progress at NASA Glenn Research Center. An early step in this program 
involves design, analysis, and testing of a Si 3 N 4 plate with internal cooling channels. In addition to the 
manufacturing and experimental work required to develop this technology, thermo-mechanical analysis of 
the Si 3 N 4 plate is needed to predict stresses that arise in the plate and to aid in the design of the plate to 
maximize cooling while minimizing resulting thermally induced stresses. Clearly, the plate geometry 
itself has limited applicability in aircraft engines, but the plate cooling channel analysis will serve as a 
template for future development of more complex shaped cooled Si 3 N 4 components. 

In addition, since functional grading of coatings (both thermal and environmental) and possibly of 
the substrate itself will be investigated in the future, the present thermal and mechanical study of the 
cooled Si 3 N 4 plate will be conducted using a recently developed higher order theory for functionally 
graded materials, referred to as HOTFGM-2D in the literature [3-7], rather than employing the traditional 
finite element analysis (FEA) technique. HOTFGM-2D offers a comprehensive approach towards 
modeling the response of material systems with different microstructural details, including certain 
advantages not available in standard displacement based finite element analysis techniques, such as: 
1) simultaneous satisfaction of displacement and traction continuity conditions between the different sub 
volumes of the spatially variable microstructure; 2) combined solution of the thermal and mechanical 
problems; 3) easy variation of cooling channel locations; 4) less mesh sensitivity; and 5) significantly less 
effort in constructing the required input. The HOTFGM-2D theory is extended herein to permit multiple 
internal boundaries of the type necessary to model a plate with internal cooling channels and now enables 
analysis of geometries with arbitrary internal and external boundaries, thus significantly expanding the 
applicability of the analysis approach to include curved cross-sections. 

The paper begins with a brief review of the previously developed higher order theory for 
functionally graded materials, including a generalization to include boundary cell classifications and 
convective thermal boundary conditions. Then two applications are addressed: the first described in 
section 3.1, is that of a flat plate with internal cooling channels and the second, described in section 3.2, is 
an axisymmetric cylinder subjected to both a thermal gradient and internal pressure. 


2.0 Higher Order Theory Review with Boundary Cell Generalization 

HOTFGM-2D is based on a geometric model of a heterogeneous composite having finite dimensions 
in the x 2 — x 3 plane and extending to infinity in the x x direction. Fig. 2. In the x 2 — x 3 plane, the 
composite is functionally graded by an arbitrary distribution of fibers of arbitrary cross section. It is 
assumed that the composite may contain one or several “windows”, which are merely holes of arbitrary 
cross section. The loading applied to the composite’s external, as well as internal, boundaries (window 
boundaries) may involve an arbitrary temperature or heat flux distribution, and mechanical effects 
represented by a combination of surface displacements and/or tractions in the x 2 — x 3 plane, and a 
uniform strain in the x x - direction which is equal to zero under plane strain conditions. 

The functionally graded microstructure in the x 2 - x 3 plane is model by discretizing the 
heterogeneous material cross-section into N q and N r generic cells in the intervals 0 < x 2 < H . 
0<x 3 < L, respectively. The generic cell (q,r) used to construct the material, highlighted in Fig. 2, 
consists of four subcells designated by the pair (/3,y), where each index /?,y take the values 1 or 2, 
which indicate the relative position of the given subcell along the x 2 and x 3 axis, respectively. The 
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indices q and r, whose ranges are q = \,2,...,N q and r = l,2,...,N r , identify the generic cell in the 

X-, — x, plane. The dimensions of the generic cell along the x 2 and -V, axes, h \ q] , fu 1 1 and l[ \ can 
vary in an arbitrary fashion such that 

H = J j (h\ q] +ltf ) ) (1) 

L=^(/; rl +/i rl ) (2) 


Given the applied thermomechanical loading, an approximate solution for the temperature and 
displacement fields is constructed based on volumetric averaging of the field equations together with the 
imposition of boundary conditions at the external and internal boundaries and continuity conditions in an 
average sense between sub volumes used to characterize the material’s microstructure. As described 
briefly in the following (and more thoroughly in references [3-7]), this is accomplished by approximating 
the temperature and the in-plane displacement fields in each subcell of a generic cell using a quadratic 

expansion in the local coordinates \ 2 P) ,\ i {Y) centered at the subcell’s mid-point. This higher-order 
representation of the temperature and displacement fields is necessary in order to capture the local effects 
created by the thermomechanical field gradients, the microstructure of the composite with internal 
windows, and the finite dimensions in the functionally graded directions. The unknown coefficients 
associated with each term in the temperature and displacement field expansions are then obtained by 
constructing systems of equations that satisfy the requirements of a standard boundary-value problem for 
the given field variable approximations. That is, the heat and equilibrium equations are satisfied in a 
volumetric sense and the thermal, heat flux, displacement, and traction continuity conditions, within a 
given cell as well as between a given cell and its adjacent neighbors, are imposed in an average sense. 

The volume discretization employed in HOTFGM-2D in the past has allowed one to represent flat 
rectangular plate geometries composed of different complex types of spatially variable material 
architectures in sufficient detail. In this work we have extended this previous capability greatly by 
considering two types of generic cells; namely, the standard internal cell (see insert in Fig. 2) and the new 
boundary cell (see Fig. 3). An important distinction between the two cells is that an internal cell is 
neighbored by other generic (internal or boundary (solid material)) cells in both x 2 and -V, directions; 
whereas a boundary cell is located on (has as a neighbor) an internal or external boundary and may or 
may not be composed of a solid material. 1 Consequently, each type of generic cell requires a different 
analysis for establishing the required equations in the unknown field variables. Also, classification of 
generic cells into either internal or boundary cells provides a convenient and computationally efficient 
means for analyzing arbitrary shaped structures with or without, internal cooling passages, cracks, or 
the like. 

Before outlining the basic analysis framework, distinction must be made between the global 
coordinates jq, x 2 , and shown in Fig. 2, and the local in-plane coordinates x,^ ,x 3 lr) used to designate 
position of each subcell (Py) within a generic cell (<:/, r) . The local and global coordinate directions 
designated by the same subscript are parallel to each other and the subcell designation ( (3y ) identifies the 
subcelfs relative position within a given cell along the global x 2 , ^coordinates. 


1 Note, when a boundary cell is not composed of a solid material, it is skipped in the analysis portion of the code. 
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2.1 Thermal Analysis 


Let the functionally graded parallelepiped shown in Fig. 2 be subjected to steady state 
temperature, heat flux or convection distributions on its external and internal surfaces. Under these 
circumstances, the heat flux, field in the material occupying the subcell (fly) of the (g,r)th generic cell 


in the region 




< tif 
- n / 




< I y r) /l must satisfy: 


dxf + dxl y> 


(fry = 1 . 2 ) 


(3) 


The components q\^ Y) of the heat flux (per unit area) vector in the subcell are derived from the 
temperature field according to the Fourier law, 

dPr) c)I 


<h 


(fir) 


= -k 


dx; 


(■) 


(i = 2,3; no summation on /) 


(4) 


where k\^ y) are the coefficients of thermal conductivity of the material in the subcell, and T m is the 
spatial temperature distribution in the subcell (fiy) of the (q, r)th generic cell, measured with respect to a 
reference temperature T refy and no summation is implied by repeated Greek letters in the above and 
henceforth. 

The temperature field in the (j3y) subcell is approximated by a second order expansion in the 
local coordinates x^\x\ y ' as follows: 


nrifir) _ 7 Afir) , —{P) r r{fty) 

* i It , 
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(5) 


where 7^^’ , which is the temperature at the center of the subcell, and 7^' 
(m y n = 0, 1, or 2 with m + n< 2) are unknown coefficients, which are determined from the conditions 
outlined subsequently. 

Given the five unknown quantities associated with each subcell (i.e., T^\ , T^ y) ) and four 

subcells within each generic cell, 20 N N r unknown quantities must be determined for a composite 
with N q and N r rows and columns of cells containing arbitrary specified materials. 

For all cells, these quantities are determined by first satisfying the heat equation (i.e., eqn (4) 
substituted into eqn (3)), as well as the first and second moment of this equation in each subcell in a 
volumetric sense in view of the employed temperature field approximation. Subsequently, for internal 
generic cells continuity of the heat flux and temperature is imposed in an average sense at the interfaces 
separating adjacent subcells, as well as neighboring generic cells. 

For boundary (internal or external) generic cells however, there are no neighboring cells in 
certain directions. Consequently, the conditions that reflect the continuity of temperature and heat flux 
between neighboring generic cells are replaced by the imposed surface boundary conditions on this type 
of boundary cells. There are three possible types of boundary conditions that can be applied currently to a 
given internal (window) or external boundary: 


1. 


The heat flux is specified at the boundary 5, namely 


4 B =( i 


applied 


( 6 ) 
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where q\ w denotes the applied heat flux at the boundary. 


2. The temperature is specified at the boundary B . namely 

71= applied 

where T\ ljed denotes the applied temperature at the boundary. 


( 7 ) 


3. Convective boundary' condition is specified at the boundary B , namely 

4=Mn ,-t.) i8 » 

where h is a parameter that varies depending upon the imposed environment (i.e., the fluid 
beyond the boundary), and denotes the environmental temperature. 


Fulfillment of these field equations and continuity conditions, together with the imposed thermal 
boundary conditions on the bounding surfaces, provides the necessary 20 N q N r unknown coefficients 

in the temperature field expansion in each (/fy) subcell. The final form of the system of 20 N q N r 
equations is symbolically represented below as: 

kT = t (9) 


where the structural thermal conductivity matrix K contains information on the geometry' and thermal 
conductivities of the individual subcells (fty) in the N q N r cells spanning the r-, and x 3 directions, the 
thermal coefficient vector T contains the unknown coefficients that describe the temperature field in each 
subcell, i.e.. 



( 10 ) 


where 

-\t T T T T (11) 

V qr ~ 1/(00)’ / (10)’ *(01)’ *(20)’ t02) J^ r 

and the thermal force vector t contains information on the boundary conditions. Solution of this eqn. (9) 
provides the unknown thermal coefficients, which, in turn give the temperature field (see eqn. (5)) in each 
subcell throughout the solid. 


2.2 Mechanical Analysis 

Given the temperature field generated by the applied inner and outer surface temperatures, heat 
fluxes, and/or convective conditions obtained in the preceding section, we proceed to determine the 
resulting displacement and stress fields. These mechanical fields are induced not only by the temperature 
distribution but also by arbitrary mechanical loading applied to the surfaces (be they internal or external) 
of the structure. 

The stress field in the subcell (/?/) of the (<?,r)th generic cell must satisfy the equilibrium 


equations 


dxi P) a*< 7) 


(7 = 2,3) 


( 12 ) 
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The components of the stress tensor o\f r \ assuming that the material occupying the subcell (/?/) of the 

(q,r)th generic cell is isotropic, are related to the strain components e^f 7 ' through the familiar Hooke's 
law: 


) S, + 2 


(13) 


where A (ft,) and are the Lame’s constants of the material filling the given subcell (/Jy), 8 {J is the 

jiffy) 

Kronecker delta and G tj is the thermal stress consisting of the products of the stiffness tensor 

components, thermal expansion coefficients, and the temperature change from a reference temperature. 
The components of the strain tensor in the individual subcells are obtained from the strain-displacement 
relations. 

The displacement field in the subcell (/3y) of the (q,r ) th generic cell is approximated by a 
second-order expansion in the local coordinates .r^ 1 and xj 7 * as follows: 



where the 40 unknown coefficients (which are the displacements at the center of the subcell) and 

0 = 2, 3) (the higher-order terms) must be determined from conditions similar to those employed 

in the thermal problem. Here, for the mechanical problem, the heat equation is replaced by the two 
equilibrium equations, and the continuity of tractions and displacements at the various interfaces replaces 
the continuity of heat fluxes and temperature. Here again we need to distinguish between internal generic 
cells and boundary cells as was done in the thermal problem. For external boundary cells, the continuity 
of displacements and tractions is replaced by the applied mechanical boundary conditions. 

Application of the above equations and conditions in a volumetric and average sense, 

respectively, produces a system of 40 N q N r algebraic equations in the unknown coefficients VL 77 | . The 
final form of this system of equations is symbolically represented by 

KU = f (16) 

where the structural stiffness matrix K contains information on the geometry and thermomechanical 
properties of the individual subcells (jSy) within the cells comprising the functionally graded material, 

the displacement coefficient vector U contains the unknown coefficients that describe the displacement 
field in each subcell, i.e„ 


where 



(17) 


«r=K»,. w mr >W 0 = 2.3) (IS) 

and the mechanical force vector f contains information on the boundary conditions and the thermal 
loading effects generated by the temperature distribution. 
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It should be noted that the above discussion pertains to the plane strain case in which the overall 
out of plane strain of the composite vanishes, i.e., £ n =0. If a generalized plane strain howevei is 
considered, £ n * 0 and is not known in advance. In such a case an additional equation must be added to 
the above system for the determination of the additional unknown out of plane strain. This equation 
expresses the fact that the overall out of plane stress must be zero, namely, O u = 0 . This stress is given 

by 


nL 4=1 ,. =1 p =i r=! 


(19) 


where } is the averaged induced stress in subcell (/Jy), which is located in the generic cell (q,r). 

The corresponding displacement in subcell (/3y) is then taken to be proportional to the out of plane 
strain, that is. 


3.0 Applications 

3.1 Internally Cooled Rectangular Plate Problem 

The primary problem that will be examined herein involves a long thin plate with ten cooling 
channels subjected to flame impingement (see Fig. 4). The substrate, which contains the cooling channels, 
consists of silicon nitride. The plate surface that is subjected to the flame has a thermal barrier coating 
(TBC) consisting of a mullite bondcoat and a porous zirconia topcoat. As the plate is free to move in the 
X] -direction (out of the plane), the problem can be treated as generalized plane strain. Similarly, since the 
plate exhibits symmetry in the x 2 - x, plane (about the x 2 -axis), we need only consider the cross-section 
geometry shown in Figs. 5 and 6. Figure 5 provides the dimensions for the baseline plate configuration 
wherein the cooling channels have identical square cross-sections and are equally spaced in the half cross- 
section. The choice of the cooling channel cross section shape was motivated primarily by ease of 
analysis at this early juncture in the cooled silicon nitride evaluation program. While in reality the 
channels will likely be circular or oval in cross-section, the square channels examined in this study are 
sufficient to evaluate the efficacy of internal cooling and to examine the effects illustrated herein. 

Figure 6 indicates the thermal and mechanical boundary conditions imposed in the HOTFGM 
analysis of the plate. Note that the boundary conditions indicated for the channel nearest to the symmetry 
boundary were applied to all cooling channels. All mechanical boundary conditions (except symmetry ) 
are traction-free. All thermal boundary conditions (except symmetry) are convective for this “baseline” 
case. The values employed for the convection coefficients ( h ) and surrounding air temperatures T m are 
admittedly not well known and were taken to be representative. For the external boundary free convection 
coefficients and the cooling channel forced convection coefficients, text book values were taken from [8] 
so as to simulate free air flow along a flat plate and forced airflow in a tube, respectively. The values 
employed for the air temperature inside the cooling channels and at the bottom and free edge ( x 3 = 0.5 
in.) boundaries are estimates provided by experimentalists [9], as is the flame thermal boundary condition 
far field temperature ( T ). The flame convection coefficient was estimated based on achieving a 
reasonable TBC surface temperature in early work on the problem [9]. On the portion of the top surface 
not subjected to flame impingement, the free convection thermal boundary condition was employed along 
with the linear T„ profile shown in Fig. 6. This profile simulates the decreasing air temperature that would 
occur as the distance from the flame increases and serves to decrease the discontinuity in the thermal 
boundary conditions. Note that a discontinuity remains at the edge of the flame as the convection 
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coefficient changes from 3.x 10 4 BTU/in 2 • s • °F to 2.04xl0' 6 BTU/in : • s • °F, see Fig. 5. Heat transfer 
due to radiation was assumed to be a higher order effect and was thus neglected. 

The temperature dependent material properties for the three plate materials (top coat, bond coat, 
and substrate) are given in Table 1. The HOTFGM analysis utilizes linear interpolation to determine the 
material properties between the given temperatures, and this analysis approach (outlined in previous 
section) involves two steps, both of which are performed during a single execution of the analysis code. 
First, the thermal problem is solved to determine the temperature field throughout the plate, immediately 
followed by the mechanical solution that determines the corresponding stress, strain, and displacement 
fields throughout the plate. 


3. 1. 1 Accuracy Assessment via Grid Refinement 

In order to determine the level of geometric refinement necessary to model accurately the 
response of the internally cooled plate, five subcell grids were constructed. The grids are shown in Fig. 7. 
The first (Fig. 7a) represents the fewest number of subcells that may be used, given the geometry of the 
plate and the top surface boundary conditions; as HOTFGM2-D requires the geometry of the material as 
well as the cooling channels to be described by generic cells, which contain 2x2 subcells. Note that the 
subcells within a given generic cell may all contain different materials if desired. This least refined grid 
consists of 8 subcells in the ^-direction and 24 subcells in the x, -direction. In Fig. 7b, the number of 
subcells in each direction has been doubled to yield a grid of 16x48 subcells. For Figs. 7c - 7e, the grid 
increases to 26x74, 48x136, and 96x272 subcells, respectively. As required for any consistent mesh 
refinement study, each refined grid contains the previous less refined (coarser) grid; that is, simply further 
subdividing the previous grid subcells forms each successive refinement. 

Model results in the form of contour plots for the temperature and stress fields in the plate are 
shown in Figs. 8 - 12 for the five different grid refinements. Note that, in the contour plots, the average 
value of the quantity (temperature or stress) for each subcell is plotted. Thus, variation of the quantities 
within the subcells is not depicted. From Fig. 8, it is clear that the 8x24 grid is not sufficiently refined to 
allow accurate prediction of the temperature field in the plate. Comparing Figs. 8a and b indicates that a 
significant change in the predicted temperature field has occurred by refining the grid from 8x24 to 1 6x48 
subcells. The maximum average subcell temperature (which occurs in the TBC) is slightly higher in the 
16x48 case, but the temperature in the bulk of the plate is significantly lower (typically more than 100 
°F). Thus it appears that the effect of the coarser 8x24 mesh is to lower the amount of cooling that occurs 
in the plate. 

Comparing Figs. 8b, c, d, and e, it is clear that further refinement does not have a significant 
effect on the temperature field for the bulk of the plate. Rather, the further refinement mainly affects the 
temperature in the TBC directly under the flame. The maximum average subcell temperature, which 
occurs in this region, increases from 2552 °F (16x48, Fig. 8b) to 2638 °F (26x74, Fig. 8c) to 2669 °F 
(48x136, Fig. 8d) to 2672 °F (96x272, Fig. 8e). This trend can be attributed to the fact that an average 
temperature is plotted for each subcell. For the less refined grids, subcells are larger, and thus the 
averaging occurs over a larger area. The subcell results are then “smeared” to a larger extent, and some 
concentrations are not depicted. As will be shown later, the point-wise temperature at the very top of the 
plate is actually highest in the case of the 8x24 grid. The minute differences in the temperature field 
between the 48x136 case and the 96x272 case indicate that, for the thermal problem, the 48x136 grid can 
be considered sufficiently refined. 

Examining the contour plots associated with the extreme grid refinements for the four present 
stress components" (CT 1]( 0 22 , <J i3 , and CJ 23 ), we see that grid refinement has a similar effect on 
mechanical behavior, see Figs. 9 - 12. In that, greater refinement does not significantly affect the general 


Note (T|2 and On are zero due to the generalized plane strain condition. 
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character of the stress component fields; rather it affects the magnitude and geometric details of the stress 
component concentrations. This can be clearly observed by comparing the 8x24 grid results (Fig. 9a - 
12a) with the 96x272 grid results (Fig. 9b - 12b). For example, comparing Figs. 9a and 9b indicates that 
the general geometric nature and magnitude of the out-of-plane (a*i direction) normal stress, <Xu, field are 
similar given either extreme of grid refinement. However, examining the concentration that occurs near 
the free edge in the region of the bondcoat/substrate interface indicates that the 8x24 grid is insufficiently 
refined to capture the fine details of the stress field. The same is true for the remaining normal stress 
components, the peel stress, (shown in Fig. 10) and the in-plane ( a\ -direction) normal stress, 
(shown in Fig. 11). In the case of the in-plane shear stress, 0\ v see Fig. 12, the concentration evident in 

12b is so localized that the 8x24 grid (illustrated in Fig. I2a) can hardly begin to capture it. However, as 
the level of refinement increases, the geometric details and magnitude of the concentration s) are better 
approximated. For all stress components only a slight difference is discemable between the results of the 
48x136 grid and those of the 96x272 grid. Thus, we conclude that, as was the case for the thermal 
problem, the 48x136 grid is sufficiently refined to model the mechanical response of the plate as well. 

Additional details regarding the effects of grid refinement are observable by plotting the field 
quantities along certain lines within the plate. Figure 13a illustrates the through thickness temperature 
(along the a\ -direction) directly under the flame (i.e., at the center of the plate along the plane of 
symmetry). Clearly, the plotted temperature field converges rapidly as the level of refinement increases. 
Note that in x-y plots such as Fig. 13a, the actual values of the field variables are plotted at five points 
within each subcell. Recall that in the contour plots, only average subcell values were plotted. This 
explains why, in the contour plot (see Fig. 8), the 8x24 grid appeared to have the lowest maximum TBC 
temperature, while in Fig. 13a, the 8x24 grid has the highest TBC temperature. Alternatively, Fig. 13b 
illustrates the horizontal (i.e., along the ^-direction) thermal profile along the top and bottom of the plate. 
Clearly, the temperature field has converged for the 48x136 and 96x272 grids. However, for the 16x48 
and 26x74 grids a perturbation is present in the temperature field along the top of the plate in the region 
of the flame edge. This perturbation is not present in the 48x136 grid results indicating that the further 
refinement eliminated the problem. 

Figure 14a shows the predicted horizontal (x 3 -direction) displacement (M3) through the thickness 
of the plate along the free edge. The plotted results show the degree to which the plate expands laterally. 
Once again it is evident that the response of the 48x136 and 96x272 grids have converged. Figure 14b 
shows the predicted vertical (x, -direction) displacement (u 2 ) along the top and bottom surfaces of the 
plate. Note that at the bottom of the plate at the symmetry boundary, a pinned boundary condition was 
employed (see Fig. 6), making this the origin for displacements. Consequently, the plotted results show 
the degree to which the plate bends. In this case all grid refinements beyond 8x24 appear to be reasonably 
well converged. 

From the contour plots of the individual stress components (Figs. 9 - 12) it is apparent that stress 
concentrations occur (as expected) in the region of the bondcoat/substrate interface near the free edge due 
to dissimilar materials. Figure 15 provides plots of the stress components along the bondcoat/substrate 
interface (in the A' 3 -direction). For the out-of-plane stress (<J U ) component, Fig. 15a shows that while the 

8x24 grid captures the general nature of the plotted CX n field, jumps occur at the subcell interfaces and the 
location of the maximum occurs further within the plate (away from the free edge) as compared to the 
more refined cases. Refinement of the grid improves these shortcomings, and the results from the 48x136 
and 96x272 grids are nearly identical. However, small jumps are still present in the 48x136 curve, and the 
maximum is slightly lower as compared to the 96x272 curve. Note that, at the free edge (*3 = 0.50 in.), 
<j u is not a traction component and thus should not vanish at this point, as illustrated in Fig. 15a. 

The peel stress, CJ 22 , along the bondcoat/substrate interface is similar to the out-of-plane normal 
stress (<J n ) in terms of the effect of grid refinement, as increased refinement leads to increased 
smoothness of the plotted curve and better representation of the concentration. Fig. 15b. Note that, at the 
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free edge (x, = 0.50 in.), the peel stress is again not a traction component and thus is not expected to 
vanish at this point. In fact, after reaching a low magnitude tensile peak just inside the free-edge, the peel 
stress becomes highly compressive, which should be beneficial, as it will help prevent delamination at the 
free edge. 

Figure 15c shows the in-plane normal stress, <7,,, along the bondcoat/substrate interface. Again, 
refinement improves the predictions, and the curves coming from the finer grids tend to converge. A large 
tensile concentration is present just inside the free edge, and in this localized region, there is a significant 
difference between the results obtained with the 48x136 grid and those obtained with the 96x272 grid. At 
the free edge ( x 3 = 0.50 in.), <X 33 is a traction component and thus should vanish at this point. Clearly, 
Fig. 15c shows that, while <7 33 is rapidly decreasing near the free edge, it does not vanish. This is merely 
due to the fact that boundary conditions are applied in an average sense to each subcell. Thus it is only the 
average of the in-plane normal stress along the subcell boundary face that should (and does) vanish at the 
free edge, not the in-plane normal stress at each point within a given subcell (i.e., here the point plotted is 
at the top of the subcell). 

Lastly, the in-plane shear stress, <J, 3 , along the bondcoat/substrate interface is plotted in Fig. 15d. 
As was the case with <7 33 , <7 23 reaches a maximum magnitude just inside the free edge, with the 
magnitude of the concentration again being somewhat lower for the 48x136 grid as compared to the 
96x272 grid. Also, like <7 33 , cr 23 is a traction component at the free edge, and thus should vanish. But, as 
before, since the results shown are along the topmost integration points within the subcells, and boundary 
conditions are applied in an average sense, the in-plane shear stress shown at the point does not vanish. 

Given the results presented in this section, it is reasonable to conclude that the 48x136 grid is 
sufficiently refined for the purposes of this study. Although this grid refinement does not completely 
capture the magnitudes of the highest stress concentrations in the plate, the vast majority of the predicted 
fields are nearly identical to those predicted using the 96x272 grid. Consequently, since the main purpose 
of this study is to demonstrate new technology and highlight certain effects of internal cooling, the 
48x136 grid has been deemed sufficient to 1) locate concentrations, 2) determine the effect of changing 
the plate configuration or boundary conditions on the concentrations, and 3) demonstrate that the model 
presented herein is effective. 1 The reason that it is desirable to minimize the size of the grid (while 
maintaining accuracy) is two fold; 1) the creation and alteration of the model input data becomes 
cumbersome as the grid becomes large and 2) execution time increases significantly with grid size. 
However, as described in the next section, thanks to a state-of-the-art linear equation solver now 
implemented within HOTFGM, execution times are no longer a primary driving factor. 


3. 1.2 Solution Optimization via Sparse Equation Solver 

In the version of HOTFGM employed in this study, the number of linear equations that must be 
solved is five times the number of subcells for the thermal problem and ten times the number of subcells 
for the mechanical problem (see Section 2). As the number of subcells and number of equations become 
large, execution times increase as well. Thus, it is desirable to employ the most efficient solution 
procedure possible. This section outlines the steps taken to improve the efficiency of HOTFGM by 
employing more efficient linear equation solvers. 

The solver originally employed [3] was a subroutine called LEQT1B that is part of the IMSL 
Fortran library and is available from Visual Numerics, Inc (http://www.vni.com). This subroutine uses 
banded storage to reduce memory requirements and uses L-U decomposition to factor the coefficient 
matrix. This solver however is known to be quite inefficient for sparse matrices; consequently, to improve 
the computational efficiency, a sparse solver was sought. This was motivated by the fact that even within 


' Though not shown, comparison to finite element solutions for certain configurations of the silicon nitride plate also 
indicated the accuracy of HOTFGM and the selected grid. 
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the bandwidth considered by LEQT1B, the coefficient matrix is typically over 95% sparse for the cases 
considered in this study. A free sparse solver subroutine known as Y12MAF was thus located. 
(http://www.netlib.org/yl2m). This solver proved to be significantly more efficient than LEQT1B, but 
was limited by the number of equations it could handle. Finally, a state of the art sparse linear equation 
solver known as UVSS, originally developed by NASA Langley and now distributed and maintained by 
SolverSoft Co. (http://www.solversoft.com) was incorporated. UVSS has proved to be quite efficient and 
able to handle a large number of equations. Sample execution times for the three linear equation solvers 
utilized are given in Table 2 for the five levels of grid refinement discussed in the previous section. These 
execution times are plotted (on a log-log scale) in Fig. 16. For the 26x74 grid, the UVSS solver was 
1,898 times faster than the LEQT1B solver and 50 times faster than the Y12MAF solver. Considering the 
required level of refinement for this study (48x136 subcell grid), UVSS clearly demonstrated its 
superiority as it was able to solve the problem in slightly more than one minute, whereas Y12MAF was 
incapable of solving the problem due to the large memory requirements. Compared to LEQT1B, which 
required almost a week to solve the problem, UVSS represents a significant (four orders of magnitude) 
improvement in efficiency. 


3.1.3 Effect of Cooling Channels and Flame Boundary Condition 

As previously stated, it is of significant technological importance to determine the effectiveness 
of internal cooling channels within a ceramic substrate. That is, do the channels provide a sufficient 
amount of cooling to justify their presence (considering the added manufacturing expense associated with 
them as well as any stress concentration they induce)? To this end, we have modeled an identical plate 
(with the baseline 48x136 subcell grid) with and without cooling channels. The subcell grid geometry is 
shown in Fig. 17a, with the corresponding temperature and /, and J 2 stress invariant profiles for each 
case being illustrated in Figs. 18a - 20a, respectively. Note that a more narrow scale (than that employed 
for the cooled plate) was used in Fig. 21 for the case of no cooling channels in order to highlight the fact 
that, although the temperature gradient is small (approximately 100 °F), one is still present. 

In concert with assessing the viability of cooling, we also wish to demonstrate the importance of 
the assumed thermal boundary conditions by including an additional case for comparison, in which the 
boundary condition employed for the flame has been altered. That is, the flame thermal boundary 
condition is taken (in Figs. 18c - 20c) to be a fixed temperature instead of the more realistic convective 
boundary conditions employed previously (baseline case, Figs. 18b - 20b). Specifically, in this case, the 
cooling channels are present and active, but now the temperature at the top of the TBC under the flame is 
prescribed to be 3507 °F. This temperature corresponds to that under the flame for the case of the 
uncooled plate with convective flame boundary condition (see Fig. 18a). 

The temperature results for the fixed temperature flame boundary condition case are shown in 
Fig. 18c. Note that an increased number of subcells were required in the region of the flame edge to 
eliminate perturbations in the temperature field similar to those depicted in Fig. 1 3b. Consequently, the 
size of the grid increased from 48x136 to 48x166 (see Fig. 17b). The corresponding thermally induced 
stress fields (i.e., I { and J 2 ) for the cooled convective and cooled fixed temperature boundary condition 

cases are shown in Figs. 19b, 20b, 19c, and 20c, respectively. 

To further highlight the differences among the three cases examined thus far (no cooling 
channels, baseline cooling channels with the convective flame boundary condition, and cooling channels 
with the fixed temperature flame boundary condition), detailed x-y plots of the temperature and 
mechanical fields (stress and displacement components) along specific lines of interest are shown in 
Figs. 22 - 24. The lines of interest were determined from examination of the various contour plots in 
Figs. 18-20. Examining Figs. 18 and 22, one immediately sees the major impact that internal cooling 
has on the temperature profile within the plate; in that internal cooling within the substrate itself typically 
lowers the temperature at a point by approximately 900 °F. Also, by adding the cooling channels 
(baseline case - convective flame boundary condition cooling), the TBC surface temperature is 
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significantly reduced. Thus, in the context of convective cooling, the presence of the cooling channels 
has, in a sense, increased the effectiveness of the TBC and created, even within the substrate, a noticeable 
though-thickness temperature drop. Alternatively, in the case of a fixed temperature flame boundary 
condition (Figs. 18c and 20b), the cooling channels still have a significant effect, but the resulting overall 
plate temperature remains much higher than in the baseline case, as one might expect, particularly when 
one considers the surface temperature of the TBC. The temperature profile associated with the fixed 
temperature flame boundary condition is similar in form to that of the baseline case, but because the TBC 
surface temperature under the flame was forced to remain at 3507 °F, the additional cooling introduced 
via the channels was unable to affect the surface temperature. Consequently, a much higher temperature 
field throughout the plate, as well as a higher gradient through the TBC layer itself, is produced. These 
points are reinforced by Fig. 22b. Note that the larger top to bottom temperature difference present for the 
fixed flame boundary condition case (as compared to the baseline case) is expected because convective 
cooling is driven by the difference between the temperature at the convective surface and T^ (see 
eqn. (8)). Hence, since the overall temperature is greater for this case, the temperature difference created 
by the convective cooling channels is larger, which leads to more cooling. Finally, considering the plate 
with no cooling channels, the small temperature drop that does occur is mainly through the TBC. 

Turning our attention to the mechanical response for the three cases considered, it should be noted 
that three major effects impact the thermomechanical response. First, the overall plate temperature (the 
highest being the case with no internal cooling and the lowest being the baseline case) drives the overall 
thermal expansion of the plate as well as the thermal expansion mismatch between the TBC and the 
substrate. The overall plate temperature also affects the temperature dependent material properties of the 
constituents. Second, the thermal gradient within the plate itself causes bending (Fig. 24b) and also sets 
up a gradient in the material properties. This thermal gradient is largest for the cooled plate with fixed 
temperature flame boundary condition and lowest for the uncooled plate. Third, the geometric asymmetry' 
of the plate, caused by the off-center cooling channel location and the presence of the TBC, causes 
additional bending. The degree of geometric asymmetry is identical for the two plate configurations with 
cooling channels but lower for the uncooled plate. All three of these influences work in concert to bring 
about the trends that are evident in the results presented. 

Examining the I x stress invariant (I { = <T n +(J 22 +<7 33 ) contour plots given in Fig. 19, one can 
immediately see the potential regions where brittle damage (crack initiation or void nucleation) could be a 
problem (i.e., regions of high tensile hydrostatic stress). Obviously, in all three cases this region is located 
within the substrate near the free edge at the bondcoat/substrate interface, with the uncooled plate having 
the highest value (73,550 psi) and the baseline (internally cooled) case having the lowest (48,987 psi). 
Also, it is clear that hydrostatic failure within the bondcoat is not an issue as the entire bondcoat is 
subjected to compression (with the maximum compressive stress occurring under the flame), as is the top 

coat and large portions of the substrate. Alternatively, examining the J 2 invariant ( J 2 = yl SijSij * where 
Sy are the deviatoric stress components) contour plots given in Fig. 20, one immediately sees regions 

where shear induce damage could potentially be a problem. From Fig. 20 we see that the bondcoat layer is 
subjected to the overall largest magnitudes of J 2 ; with the maximum occurring within the bondcoat near 
the centerline of the plate (under the flame). Similarly the maximum J 2 within the substrate appears to be 
underneath the flame. Consequently, given this stress overview, we have displayed the individual stress 
components as a function of horizontal (x 3 ) position within the substrate along the bondcoat/substrate 
interface, see Fig. 23, for all three cases. 

From Fig. 23a, it is clear that the largest tensile out-of-plane stress (C7 n ) occurs just inside the 
plate near the free edge. Note that the (T n concentration has the greatest magnitude for the fixed 
temperature flame boundary condition and the smallest magnitude for the plate with no cooling channels. 
Thus it appears that this tensile concentration is brought about mainly by the thermal gradient (which is 
largest for the fixed temperature flame boundary condition) rather than the overall plate temperature 
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(which is highest for the plate with no cooling channels). Alternatively, as suggested previously by 
Fig. 19a, but not specifically shown, the maximum compressive out-of-plane stress occurs within the 
bondcoat and is largest for the plate with no cooling and smallest in the cooled convective (baseline) case. 
This suggests that the overall plate temperature drives the maximum compressive <7,, stress. 

Figures 23b and contour plots (not shown) of the peel stress (CT 22 ) component illustrate that for 
all three cases considered, the peel stress is close to zero throughout the plate except in the vicinity of the 
free edge, as one might expect. Further, as noted previously, the peel stress at the free edge, located at the 
substrate/bondcoat interface is highly compressive, which may help prevent delamination. The magnitude 
of the compressive peel stress concentration is greatest for the uncooled plate and smallest for the baseline 
case and thus appears to be driven by the overall plate temperature. The maximum tensile peel stress at 
the substrate/bondcoat interface (see Fig. 23b) is also greatest for the uncooled plate. However, although 
not shown, this local maximum is not the maximum in the entire plate for any of the three cases. For in 
both cases with cooling channels, the absolute maximum occurs at the upper right comer of the cooling 
channel closest to the free edge, i.e., point A of Fig. 19c. This concentration is significantly higher for the 
fixed temperature flame boundary condition as compared to the baseline (convective) case. Alternatively, 
for the uncooled case, the maximum occurs within the bondcoat at the free edge. Note that the peel stress 
is a traction component at the substrate/bondcoat interface and is thus continuous (in an average sense) 
across this interface. 

Now considering the in-plane normal stress (<7 33 ) (shown in Fig. 23c) it is clear that, as was the 
case with the out-of-plane stress (0^,), (7 33 attains a maximum tensile magnitude in the substrate near the 
free edge along the substrate/bondcoat interface. Alternatively, as suggested by Fig. 19a but not explicitly 
shown, the maximum compressive tJ 33 occurs in the bondcoat in the region under the flame. The trend in 
the magnitudes of both of these stress concentrations appears to be driven by the absolute temperature, as 
the magnitudes are largest for the plate with no cooling and lowest for the baseline (convective) case. The 
reduction in the maximum tensile in-plane normal stress (<T 33 ) associated with the introduction of the 
cooling channels (shown in Fig. 23c) may be of significant importance as this stress component was 
reduced from 54 ksi (for the uncooled case) to 30 ksi (for the baseline case). Particularly, since the tensile 
strength of the silicon nitride substrate (at 2552 °F) is 58 ksi [10] the predicted temperature (for the 
uncooled plate) in this region of the plate is 3418 °F. Alternatively, in the case of cooling the predicted 
temperature is reduced to 1833°F and consequently the strength should be greater as well. 

In-plane shear stress (C7 23 ) results for the three cases considered thus far are illustrated in 
Fig. 23d, where, similar to the peel stress component, the shear stress is low throughout the plate except 
near the free edge where a concentration develops at the substrate/bondcoat interface. The magnitude of 
the concentration is largest for the uncooled plate and smallest for the cooled baseline (convective) case 
and thus again appears to be driven by the plate’s overall temperature. Clearly this large interfacial shear 
stress represents a potential failure mechanism for the plate that could lead to interfacial delamination, 
depending on the interfacial bond strength of the substrate/bondcoat interface. 

Finally, considering both the horizontal displacement (w 3 ) at the free edge (Fig. 24a) and the 

vertical displacement ( u 2 ) along the top and bottom of the plate (Fig. 24b) for all three cases, one can see 
the manifestation of the influence of temperature (overall magnitude, gradient or profile), geometry (layer 
configuration and thickness, channel presence and location), and material properties (temperature 
dependence, mismatch). Given the temperature profiles (overall magnitude and through-thickness 
gradient) displayed in Fig. 22, one can explain the ordering of the resulting displacement profiles given in 
Fig. 24. To do this it is important to remember the following key factors: 

1. Given the current layered system’s asymmetry (wherein the topcoat and bondcoat have higher 
CTE than the substrate, see Table 1) downward bending would occur even under a uniform 
temperature field. 
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2. This bending would increase in proportion to the overall temperature rise. 

3. Superposition of a positive through-thickness temperature gradient (i.e., higher temperature in 
the topcoat, lower temperature in the substrate) would increase bending; whereas a negative 
through-thickness temperature gradient (lower temperature in the topcoat, higher temperature 
in the substrate) would decrease bending. This effect is, however, muted by the temperature 
dependence of the Si^N 4 elastic modulus (see Table 1). which causes cooler regions to be 
stiffer and thus offer greater resistance to bending. 

4. Introduction of channels within the substrate effectively reduces the substrate’s structural 
stiffness (resistance to bending) and thus increases the amount of overall bending. 

5. The vertical location of the channels can either increase or decrease the overall system on 
bending through the effect of plate asymmetry. 

Clearly, the ordering in Fig. 24a is determined by the overall (or average) plate temperature. That 
is, the plate with no cooling channels is hottest and thus exhibits the greatest amount of lateral expansion 
(« 3 ), followed by the cooled fixed temperature flame boundary condition case, and finally by the baseline 
case, which is coolest. However, the bending profiles (measured by the displacement of the bottom of the 
plate) in Fig. 24b do not follow this simple trend. That is, based on overall temperature, one would expect 
the ordering (from least to most bending) to be: baseline, fixed temperature boundary condition, no 
cooling channels. Based on temperature gradient, this trend would be: no cooling channels, baseline, fixed 
temperature case and lastly, based solely on geometric considerations, the baseline and fixed temperature 
cases would be expected to bend more than the uncooled plate. Apparently, the observed trend in bending 
in Fig 24b is: baseline, no cooling, followed by fixed temperature case. 

Given the aforementioned facts, one can explain the observed ordering of the bending results in 
Fig. 24b in terms of the interactive nature of the three influencing factors. First let us examine the 
uncooled case, which has almost no through-thickness temperature gradient, but the highest overall 
temperature. This case exhibited more bending than the cooled (convective) baseline case that had a 
relatively significant through-thickness temperature gradient. Clearly, the effect of the higher overall 
temperature in the uncooled plate has dominated the influence of the greater thermal gradient and lower 
effective substrate stiffness in the baseline case. Next considering similar geometries, i.e., the two cases 
with cooling channels, one sees the anticipated ordering wherein the case with the larger through- 
thickness thermal gradient and overall temperature (fixed temperature flame boundary condition) bends 
significantly more than the baseline convective flame boundary condition case with its smaller gradient 
and overall temperature. Finally, comparing the fixed temperature flame boundary condition case (which 
exhibits the most bending) we see that here the effects of the higher thermal gradient coupled with those 
of the reduced effective substrate stiffness dominate over the influence of the higher overall temperature 
in the uncooled plate. These results clearly demonstrate the important interaction of temperature, 
geometry, and material properties and point to the positive impact of reducing the overall temperature 
through cooling on the plate’s deformation behavior. 

In summary, this section clearly illustrates the benefits of introducing cooling channels into the 
plate. The temperature throughout the plate can be significantly reduced, which, as discussed previously , 
can significantly improve the fatigue life of the plate (or a similar component). But the benefits go further. 
The presence of the cooling channels also reduces the magnitudes of the stress components (except for the 
tensile <T U concentration) and displacements. These effects can also be expected to improve the failure 
and fatigue life characteristics of a cooled plate or similar component. Furthermore, from a modeling 
viewpoint, the importance of employing realistic thermal boundary conditions (particularly when 
evaluating the effectiveness of internal cooling) was demonstrated. In that, if a prescribed surface 
temperature instead of convective boundary conditions is used to model the effect of the flame on the 
plate, the impact of cooling channels on reducing surface temperatures under the flame is artificially 
reduced. Consequently, an unrealistically high temperature field throughout the plate, which also leads to 
greater predicted stresses and displacements, results. As for the effects that contribute to the mechanical 
behavior of the plate (overall temperature, temperature gradient, geometric asymmetry, and material 
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properties), the results in this section indicate that the overall temperature is usually the most significant 
of these factors, with the temperature gradient followed by the degree of plate asymmetry being less 
important. 


3.1.4 Effect of Cooling Channel Location 

Here, the effect of moving the cooling channel locations with respect to each other and the plate 
boundaries is examined. HOTFGM is ideally suited for conducting such parametric studies as 1 ) the input 
is considerably more simplistic than that of similar finite element models, and 2) the thermal and 
mechanical analyses are performed sequentially within a single execution of the code. Hence, employing 
the HOTFGM analysis code significantly reduces turnaround time for a large number of cases. When 
defining the allowable cooling channel locations, a geometric constraint was imposed so as to prevent the 
channels from becoming too close to a plate boundary or each other, thereby avoiding excessive thermal 
or mechanical gradients from being introduced. A minimum distance of 0.01 inches was used for this 
purpose. 

Horizontal Channel Location 

The influence of perturbing the channel positions along the same horizontal line is examined first. 
Four cases are considered: 1) the baseline, 2) uniformly distributed channels but shifted left, 3) uniformly 
distributed channels but shifted right, and 4) channel spacing increasing linearly from left to right. 
Thermal and mechanical contour plots for these cases are provided for comparison in Figs. 25 - 27. 
Recall that the left boundary (center line) of the plate has symmetric boundary conditions, thus the 
placement of the left-most cooling channel a distance of 0.005 inches from this boundary actually 
represents a channel spacing of 0.01 inches. Figures 28 — 30 provide x-y plots and a bar chart for the four 
plate configurations presented in this section. 

Comparing the plate temperature field (Figs. 25 and 28) among the four cases, it is clear that the 
horizontal location of the channels has a noticeable effect. By moving the channels to the left (Fig. 25b), 
the temperature in the substrate under the flame is reduced. On the bottom surface, at the midpoint of the 
plate, the temperature is reduced by nearly 100 °F compared to the baseline case (Fig. 25a). Further, the 
temperature on the plate surface (surface of the topcoat) directly under the flame increases slightly as a 
result of shifting the channels to the left, as does the temperature at the bondcoat/substrate interface 
(Fig. 28c). At the free edge of the plate, the temperature is higher than the baseline case by approximately 
18 °F when the channels are shifted to the left. 

Alternatively, if the channels are shifted to the right (Fig. 25c), the temperature in the plate near the flame 
increases significantly compared to the baseline case (from 2058 °F to 2123 °F at the bottom surface of 
the plate directly under the flame), while near the free edge, the temperature decreases slightly. Thus, the 
results of shifting the channels horizontally indicate that by moving the channels to the left, which is 
closer to the heat input, the cooling of the plate is improved. Conversely, moving the channels to the right 
(away from the heat source) is detrimental to the overall cooling of the plate. The exception is near the 
free edge, where the stress concentrations are highest. 

To capitalize on the improved overall cooling provided by the horizontal shift toward the heat 
source (left), yet diminishing the temperature increase near the free edge, one could employ a linear 
channel distribution (Fig. 25d). By placing the first two channels close to the heat source, the overall 
cooling of the plate is improved compared to the baseline case (see Fig. 28). However, by allowing the 
right-most channel to be positioned closer to the free edge as compared to the uniformly distributed left- 
shift case (Fig. 25b), the temperature near the free edge is increased only slightly (approximately 1 .5 °F) 
as compared to the baseline case. 

Turning our attention to the resulting thermally induced invariant stress profiles (i.e., /, and J 2 
depicted in Figs. 26 and 27, respectively) one immediately sees how the regions of high hydrostatic 
tensile stress move with channel placement and how the overall maximum compressive and tensile states 
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are modified. Note however that the sign of the stress state within a given layer is not fundamentally 
changed (due to channel placement) from that of the baseline discussed previously. Furthermore, we see 
that the magnitude of the hydrostatic stress, /,, is increased and shifted by moving the channels to the 
right and decreased by moving the channels to the left. Also, the presence of the cooling channels near the 
free edge is clearly detrimental in that both /, and J 2 stress invariants as well as the (T n component 
(Fig. 29a) are increased and a large damaging stress concentration arises in the lower half of the plate at 
the free edge. Similar concentrations also arise internally in regions of the substrate between the channels. 
Finally, neither the in-plane normal stress, <7 33 , (Fig. 29c) nor the in-plane shear stress, 0\ 3 , (Fig. 29d) 
are significantly affected by the horizontal placement of the cooling channels. 

In Fig. 30 the displacement profile for both the free edge and the top and bottom surface of the 
plate are illustrated. As mentioned previously, depending on the particular application of a given 
internally cooled component, the amount of displacement may or may not be of criticality to achieve 
design performance. With this in mind. Fig. 30 illustrates that the plate with the channels shifted to the 
right would perform more poorly, as it exhibits the greatest amount of extension and bending. Both the 
plate with the channels shifted left and the plate with the linear channel spacing exhibit decreased 
extension in the .v 3 -direction yet increased bending compared to the baseline case. Among these three 
configurations, a trade-off is present between the beneficial effects of the increased cooling and reduced 
extension and the detrimental effects of the increased thermal gradient and bending. Shifting the channels 
to the right produces only detrimental effects and suggests that placing cooling channels away from the 
heat source (near the plate’s free edge) should be avoided. 

Vertical Channel Location 

Next, the influence of altering the vertical position of the cooling channels is examined. As in the 
previous horizontal location study, contour plots for the comparison of the resulting temperature and 
stress profiles are shown in Figs. 31 - 33. The positioning of the channels within the plate is most easily 
seen by examining the temperature profiles given in Fig. 3 1 wherein again four cases are examined: 1 ) the 
baseline (Fig. 31a), 2) the case in which the channels have been moved up toward the TBC layer as far as 
possible while still maintaining the 0.01 inch minimum distance from the bottom of the TBC layer as 
discussed previously (Fig. 31b), 3) a configuration in which the channels have been moved as far 
downward as allowable (Fig. 31c) and 4) a configuration which produced the minimum bending 
discussed later in this section (shown in Fig. 3 Id). Obviously, altering the vertical position of the cooling 
channels affects the temperature field within the plate in a similar (yet less significant) manner then did 
shifting the horizontal channel position. However, the vertical movement of the channels has the 
additional effect of changing the through-thickness asymmetry of the plate and thus influencing the 
bending of the plate as well through this additional asymmetry mechanism. Consequently, as will be 
shown later in this section, it is possible to determine the vertical location of the channels such that the 
amount of plate bending is minimized. 

Again from Fig. 31, it is clear that moving the channels up toward the flame does lower the 
overall plate temperature somewhat even though the topcoat temperature increases. For example, at the 
midpoint of the lower plate surface (beneath the center of the flame), the temperature is reduced by 8.7 °F 
compared to the baseline case; whereas along the lower surface at the free edge, this reduction is 21.4 °F. 
Alternatively, by moving the channels down (Fig. 31c), the resulting temperature at these locations is 
increased by 15.7 °F and 14.0 °F, respectively, as compared to the baseline configuration. Perhaps most 
interesting is the local increase in temperature that occurs directly under the flame when the channels are 
moved up (see Fig. 34c). The temperature increase within this region is a result of the fact that the volume 
of the substrate material into which the heat from the flame can flow has been reduced. However, as 
evidenced by the lower plate temperature outside this region, heat flow into the cooling channels 
increases as well. Thus, if the material can withstand the higher local temperature directly under the 
flame, the configuration shown in Fig. 31b could be beneficial. 
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Examining the contours of hydrostatic (/,) and von Mises ( J 2 ) stress in Figs. 32 and 33, one 
immediately observes that moving the cooling channels upward (toward the heat source), see Figs. 32b 
and 33b, results in a significant increase in both tensile 7, and J 2 stress invariants within the substrate 
material just below the bondcoat. Moving the channels downward (Figs. 32c and 33c) lowers both of 
these invariant quantities. The maximum tensile hydrostatic state appears always to be in the vicinity of 
the free edge at the substrate/bondcoat interface. The maximum J 2 invariant is always in the bondcoat 
layer under the flame; whereas the maximum J 2 within the substrate typically occurs near the free edge at 
the substrate/bondcoat interface. The only exception is when the channels are moved upward. 

Figure 35a indicates that the out-of-plane normal stress (Oj,) component concentrations are 
somewhat increased by moving the cooling channels vertically upward, whereas moving the channels 
downward has very little effect on this component. Alternatively, the tensile peel stress. <7,,, 
concentration near the free edge is unaffected by vertical channel movement, although the beneficial 
compressive magnitude at the free edge is significantly increased by moving the channels up (see Fig. 
35b). Considering the in-plane normal stress, <J 33 , distribution it is evident that this component is only 
slightly affected by moving the channels down, whereas moving the channels up significantly increases 
the magnitude of this stress component (see Fig. 35c). Similarly, moving the cooling channels up 
increases the magnitude of the in-plane shear stress. <j ;3 , concentration near the free edge at the 
bondcoat/substrate interface; whereas moving the channels downward has only a small effect on this <7 23 
component (see Fig. 35d). 

Figure 36 shows the displacement results for the cases in which the vertical location of the 
channels has been altered while maintaining the baseline horizontal positions. As Fig. 36b indicates, 
moving the channels both up and down results in increased bending of the plate. The fact that both 
configurations produce more bending than the baseline case seems counterintuitive since from a purely 
mechanical standpoint resistance to downward bending should be increased by moving the channels up, 
whereas moving the channels downward should decrease the resistance of the plate to bending. However, 
one must realize that two competing effects are at work; 1) the through-thickness temperature gradient 
and 2) the plate asymmetry (due to vertical channel location). Additionally, it must be remembered that 
the temperature gradient is significantly impacted by the cooling channel placement (asymmetry) in that if 
the channels are moved upward the through-thickness gradient increases whereas if the channels are 
moved toward the bottom of the plate the gradient is decreased. Note that the sign of the gradient is 
always maintained, thus increasing (or decreasing) the gradient only magnifies (or lessens) the downward 
bending of the plate. Results indicate that the thermal gradient is the more dominant driving force for 
bending as compared to the plate asymmetry. Consequently, it should be possible to find a vertical 
channel location that minimizes the induced bending. Figure 36c shows the displacement of the bottom 
plate surface for five different vertical channel locations. If the deflection of the lower right comer of the 
plate is used to quantify the amount of plate bending, we see in Fig. 36d that, indeed, a minimum bending 
configuration exists when the bottom of the cooling channels are located 0.022 inches from the bottom of 
the plate. Note that this configuration fortuitously is only slightly different from the baseline case and 
consequently produces only a slight reduction in the amount of bending as compared to the baseline case. 
This minimum bending study is intended to be an illustrative example, rather than an important design 
conclusion for the particular internally cooled plate investigated in this study, as optimization considered 
only vertical channel placement. The influence of combining vertical and horizontal channel placement 
will be studied in the next section. 

Combined Horizontal and Vertical Channel Location 

As a final set of illustrations of channel location effects, we examine two cooling channel 
configurations in which the channels have been shifted both horizontally and vertically. Fig. 37b shows 
the first configuration in which the channels are shifted as far left and as far up as possible considering the 
self-imposed 0.01 inch minimum spacing. This configuration represents an attempt to maximize the 
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overall cooling of the plate. Alternatively, Fig. 37c shows a configuration in which the left-most channel 
has been placed as far up and to the left as possible, with the remaining channels having their horizontal 
positions linearly increased in the .redirection. Note, the second channel (from the left) is placed 
(vertically) in the center of the substrate, while the remaining three channels’ vertical positions have been 
moved towards the bottom of the plate. This configuration represents an attempt to capitalize on the 
increased cooling associated with placing a cooling channel close to the flame-impinged surface while 
diminishing some of the negative effects associated with the previous configuration. 

As shown in Figs. 37 and 38, the two new channel configurations offer noticeable temperature 
reductions throughout the plate with the exception of the region above the first channel directly under the 
flame (the reason behind this local heating was discussed previously). In the first configuration (Fig. 37b), 
the temperature at the midpoint of the bottom of the plate is reduced by more than 1 00 °F compared to the 
baseline case (Fig. 37a). Similarly, as was the case when the channels were just moved left (Figs. 26b and 
29a), the magnitude of the out-of-plane normal stress, <7 n , concentration near the free edge is also 
reduced for this configuration (Fig. 41a); even though the magnitude of I x (see Fig. 39) and specifically 
<J H is increased elsewhere in the plate, particularly above the channels and at the free edge. In contrast, 
for the linearly spaced and vertically arranged channels, the tensile hydrostatic stress maximum is reduced 
slightly as compared to the baseline case, however now a tensile concentrations arises for this 
configuration near the top left comer of the second cooling channel from the left. This concentration is 
potentially problematic as this is also an area of high temperature (lower material strength). 

The first significant detrimental effect associated with moving the channels up and to the left is in 
the in-plane normal stress (cr 33 ) field, as the magnitude of the highest <J 33 concentration (located at the 
bondcoat/substrate interface near the free edge) is increased noticeably compared to the baseline case 
(see Fig. 41c). Also, the magnitude of this stress component is increased throughout the region of the 
substrate located above the cooling channels thereby potentially increasing the probability of failure due 
to the larger volumetric sampling of flaws. The alternative configuration with linear spacing and vertical 
arrangement reduces the in-plane normal stress near the free edge (as compared to the previous case), but 
now a concentration arises at the upper right comer of the left-most channel (not shown). Note however, 
that this new concentration may potentially be more detrimental to the plate than the higher concentration 
in the plate with the previous arrangement as this new concentration is in a region of higher temperature 
and thus lower strength. 

The in-plane shear stress, <J 23 , in the plate behaves in a similar manner to that of the <7 33 

component in that the magnitude of stress is increased by moving the channels up and to the left as 
compared to the baseline configuration. Whereas, for the linear horizontal and vertical channel 
arrangement, the magnitude decreases somewhat. Once again, a concentration arises at the upper right 
comer of the left-most cooling channel for this configuration (not shown), but in the case of 0\ v the 
magnitude of this new concentration is not very large. Note that along the bondcoat/substrate interface 
(Fig. 4 Id) no real difference is observed in the <J 23 results. 

Figure 42 shows the most significant consequence of locating the cooling channels as far up and 
to left as possible, in that both displacement components are increased significantly as compared to the 
baseline case. The amount of bending undergone by the plate is almost doubled compared to the baseline 
case. As discussed previously, this may or may not be of concern for a particular application. However, 
recall that if the plate were fixed rather than free and thus unable to bend or extend, the stresses in the 
plate would be increased, perhaps significantly, due to this channel configuration. As Fig. 42 illustrates, a 
great deal of the increased bending associated with moving the channels up and left can be eliminated by 
employing the linearly spaced and vertically arranged cooling channel configuration. Clearly to obtain the 
best overall design a full thermal-mechanical optimization study (taking into account an appropriate 
failure criterion) must be undertaken due to the inherent complex thermomechanical coupling present in 
the problem. 
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3.2 Modeling Curved Cross-Sections 

Now with the introduction of boundary cells, the modeling capability of the higher order 
micromechanics theory, HOTFGM, has been extended to include arbitrary-shaped structures as well as 
those containing internal cooling passages. To illustrate this new capability, the thermo-elastic behavior 
of components with curved cross-sections will be examined. First, the thermo-elastic behavior of a thick- 
walled cylinder subjected to internal and external pressure and thermal convection loading will be 
modeled. In this way, the HOTFGM results can be validated against the corresponding analytical solution 
results derived in the Appendix. 

For consistency a cylinder composed of silicon nitride was examined. The temperature 
independent material properties are given in Table 1 at 77 °F and the remaining required analysis input 
parameters describing the problem are given in Table 3. The geometry and boundary conditions 
associated with the cylinder problem are shown in Fig. 43a, while Fig 43b show the differential volume 
element employed in the analytical solution. The HOTFGM geometric idealization of the cylinder being 
analyzed is shown in Fig. 44. Due to symmetry considerations, only one quarter of the cylindrical cross- 
section was actually modeled and analyzed. As Fig. 44 indicates, curved boundaries must be represented 
with steps, as the current HOTFGM analysis methodology still requires rectangular subcells. The impact 
of this geometric constraint on the accuracy of the analysis will be discussed and illustrated later. 

Accurately reproducing the required cylindrical coordinate boundary conditions in the inherently 
Cartesian HOTFGM framework is challenging but feasible. Simply applying a normal stress on every 
boundary face in Fig. 44b is not sufficient, as the actual pressure boundary conditions require a radial 
stress to be applied, see Fig 43a. Consequently, for use as boundary conditions in HOTFGM, the in-plane 
stresses in Cartesian coordinates must be related to the in-plane stresses in cylindrical coordinates through 
the known transformation equations: 

a 2 = a r sin 2 (0) + <J 0 cos : (f?) + 2 a r6 sin(0)cos(0) 

(7 3 = o r co$ 2 (6) + <7 e sin 2 (0)- 2<7,0 sin(0)cos(0) (22) 

( 7 23 = ((J r ~a^)sin( 0 )cos( 0 )-f-( 7 r0 (cos 2 ( 0 )-sin 2 ( 0 )) 

On the cylinder boundaries, c 7 r is the only stress component present, but as eqns. (22) indicate, in general 
this leads to the imposition of all three Cartesian stress components. In HOTFGM, however, it is not 
possible to apply all three components on each face of the boundary since only one of the in-plane normal 
stress components represents a traction on each face. Toward this end, eqn. (22) was used to determine 
the Cartesian stress components for each boundary subcell face, given the value of <7 r and the value of 6 
for the midpoint of the boundary subcell face. The shear stress and one of the two normal stresses were 
then applied as the boundary conditions to simulate the internal and external pressure. 

A similar procedure was employed relative to thermal boundary conditions; since applying 
convection normal to each boundary face in HOTFGM is unrealistic, as the actual heat flux due to 
convection occurs in the radial direction, not in a given Cartesian direction. However, unlike stress, heat 
flux is a first-order tensor and thus must obey the applicable transformation equations: 

q, -q r sin(0) + <7, cos(0) 

(23) 

q^=q r cos(d) + q e sin(0) 

For the present cylindrical problem, only q r is present on the boundary and consequently is applied (in 
the analytical solution) through convection as shown in eqn (A.8) of the appendix. As this equation 
shows, for convection, the heat flux q is proportional to the convection coefficient h. Thus, we can 
replace eqns. (23) with. 
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( 24 ) 


h 2 = h r sin(0) 
h 3 = h r cos(0) 

and use these equations to determine the appropriate value of the convection coefficient to apply to each 
boundary subcell face given the direction of the normal to that face (2-direction or 3-direction); that is, the 
actual value of the radial convection coefficient for the boundary and the value of 0 for the midpoint of 
the subcell face. 

The actual subcell grid geometry for the HOTFGM idealization of the cylindrical problem 
analyzed is shown in of Fig. 44b. The grid consists of 100 subcells in each direction, but since only 2828 
of the total 10,000 subcells are material occupied, the execution of the HOTFGM code is still very rapid. 
Note that the boundaries located on and along the x 2 and a' 3 axes have symmetric boundary conditions 
imposed, thus the entire cylinder is, in effect, being modeled. Figure 45 shows the temperature field for 
the cylinder predicted with HOTFGM. Recall that the simulated thermal problem involves a cool internal 
fluid with a high convection coefficient and a hot external fluid with a lower convection coefficient (see 
Table 3 and Fig. 43a). Thus, the temperature profile of the cylinder must increase from the inner to outer 
surface. The temperature field shown in Fig. 45 also shows little angular (6) dependence (i.e., the 
temperature remains essentially constant at a particular radial position as the angle changes). Recall that, 
due to the axisymmetry of the problem, the analytical thermal solution depends only on radial position 
(see eqn (A. 7) of the Appendix). 

The analytical and HOTFGM solutions are compared in Fig. 46, where the temperature versus 
radial position is shown given an angle of approximately 84.5 degrees form horizontal, see Fig. 44b. The 
value of the radial position of every point along this line was determined from the Cartesian position of 
the subcell centroid via rotation equations. Note, this line is a sufficient distance from the left boundary 
(Xy = 0) to avoid any local effects (disturbances) that this boundary might cause. It is clear from Fig. 46 
that agreement between HOTFGM and the analytical thermal solution is excellent. 

Figure 47 shows individual stress component contour plots resulting from the HOTFGM 
thermomechanical problem, involving both the applied temperature field and the high internal pressure 
and lower external pressure (see Table 3 and Fig. 43a). The cylindrical coordinate stress components 
plotted in Figs. 47b - d were determined via rotation equations, given the Cartesian coordinates of the 
subcell centroids and the values of the Cartesian coordinate stresses. As was the case with the temperature 
field, the symmetry of the problem dictates that the resulting stresses should depend only on radial 
position, not on angular location. The longitudinal (along the long axis of the cylinder) stress field 
(Fig. 47a) exhibits this angular independence reasonably well, although along the inner surface of the 
cylinder, the longitudinal stress does not remain completely constant. However, this discrepancy is 
expected since the mechanical boundary conditions imposed to simulate radial pressure were not 
absolutely correct as only one of the two normal stress components needed to represent perfectly the 
radial pressure could be applied. 

Figure 47b shows the radial stress component in the cylinder as predicted by HOTFGM. Note that 
the magnitude of this stress is smaller than that of the longitudinal stress, thus smaller variations appear 
more significant. The basic character, however, of the lack of angular independence is present in the 
radial stress field given the inaccuracies associated with the application of the pressure boundary 
conditions being present once again. Similarly, it also appears that the presence of the stair step (subcell 
comers) along both inner and outer surfaces of the cylinder leads to perturbations in the stress field. The 
hoop stress (Fig. 47c) appears to be better behaved. Note that the magnitude of the hoop stress is larger 
than both that of the radial and axial stress, and the inaccuracy at the inner cylinder surface persists. 
Figure 47d shows (in cylindrical coordinates) the in-plane shear stress field. Once again, from symmetry 
arguments, the in-plane shear strain, and thus the in-plane shear stress, must be zero. Clearly, from 
Fig.47d, this shear stress component is not zero, although it is relatively small throughout most of the 
cylinder. As before, some deviation occurs at the inner and outer surfaces of the cylinder. 
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As in the pure thermal case, in Fig. 48 the HOTFGM stress results are compared with the 
corresponding analytical solution along a given angle (i.e., 84.5 degrees) from the horizontal axis. The 
overall agreement between the analytical solution and HOTFGM once again is excellent, even though a 
small deviation does occur in the hoop stress near the inner surface. Note that the analytical solution shear 
stress plotted in Fig. 48 is identically zero. Clearly, the shear stress coming from HOTFGM is quite close 
to zero, except near the inner surface. Again, this deviation is expected given the method employed to 
apply the radial pressure boundary conditions. Note that the deviation of the HOTFGM solution from the 
analytical solution near the inner surface could increase if the results were plotted along a line that passed 
through one of the larger perturbations shown in Fig. 47. 

It is clear from the above cylindrical analysis example that HOTFGM (although originally 
developed to model flat plates) with the incorporation of boundary cells can now be used to accurately 
model components with curved surfaces. Thermal problem solutions for such structures can be considered 
accurate, even in the vicinity of the curved boundary, whereas mechanical problem solutions can be 
considered only accurate in general, where in the immediate vicinity of the curved surfaces, some 
perturbations should be expected. 

As a final illustrative example, we will consider a technologically important thermal problem 
with curved surfaces that does not have an exact analytical solution: a hollow airfoil cross-section. The 
geometry is shown in Fig. 49a, where only the upper half of the airfoil need be analyzed due to the 
employment of symmetric boundary conditions. The actual subcell grid employed consists of 30x154 
subcells and can only be used to perform a thermal analysis as the subcell grid is not sufficiently refined 
to accurately model the mechanical response of the airfoil. The thermal boundary conditions imposed 
consist of a convective flame along the flat leading edge boundary with free convection on all other 
external boundaries. Internally, we consider two convective boundary conditions: free convection to 
model a hollow airfoil with no internal cooling and forced convection to model a hollow airfoil with 
internal forced air cooling. The values of the convection coefficient and employed for the flame are 
those used previously, that is, 0.0003 BTU/in. : • s °F and 3600 °F, respectively. For the free convective 
boundary conditions, these values are 2.04xl0’ h BTU/in.“ • s ■ °F and 1292 °F, and for the forced 
convective boundary conditions these values are 3.87xl0' 5 BTU/in." ■ s • °F and 1292 °F. To account for 
the surface curvature along the stepped boundaries, the convection coefficient values were multiplied by 
the sine or cosine of the angle formed by the step, depending on whether the face was normal to the x 2 or 

V; direction. 

The resulting temperature profile for the case with no internal cooling is shown in Fig. 49b, while 
the case with internal cooling is shown in Fig. 49c. Clearly, internal cooling has a significant effect, in 
that the temperature at the airfoil leading edge is reduced from 3228 °F to 2450 °F and that at the trailing 
edge is reduced from 2874 °F to 1614 °F. Consequently, armed with the current HOTFGM capabilities 
(i.e., generic boundary cell formulation and convective boundary condition capability) there is no reason 
that an airfoil with inter-wall cooling channels, like those shown in Fig. 50, could not be modeled as well. 
Such a study will be the topic of future study as experimental results become available. 


4.0 Conclusion 

This paper presents a brief summary of a Cartesian-based higher order theory for functionally 
graded materials which has been extended to include boundary (both internal and external) cells that 
enables the thermoelastic analysis of arbitrary shaped, actively cooled, structures with spatially varying 
microstructures in two orthogonal directions. A new convective thermal boundary condition and sparse 
solver were also added. Both the flat plate and cylinder applications discussed herein were only composed 
of monolithic rather than functionally graded layers of materials and therefore did not exercise the full 
capability of the theory. In particular, the viability of introducing cooling channels into a ceramic plate 
with a thermal barrier coating was examined in detail. It was demonstrated that the temperature 
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throughout the plate can be significantly reduced by the introduction of internal cooling channels, thus 
significantly improving the fatigue resistance of the plate (or a similar component). But the benefits go 
further, in that the presence and location of the cooling channels also reduce the magnitudes of the stress 
components (except for the tensile On concentration) and displacement field, which in turn can positively 
impact the failure and fatigue life characteristics of a cooled plate or similar component. As for the effects 
that were shown to contribute to the mechanical behavior (deformation) of the plate (overall temperature, 
temperature gradient, geometric asymmetry, and material properties), the results indicate that the overall 
temperature is typically the most significant of all factors, with the temperature gradient followed by the 
degree of plate asymmetry as the next most important factors. Specific observations regarding the 
influence of cooling hole placement are as follows: 

• Channel placement has approximately an order of magnitude less impact on temperature 
reduction as compared to the existence or nonexistence of cooling channels. 

• Moving cooling channels horizontally toward (away from) heat source improves (is 
detrimental to) the overall cooling of the plate, relative to uniform cooling channel spacing; 
whereas extension and bending are increased. Grading the horizontal spacing of the cooling 
channels will increase overall cooling while reducing extension and other detrimental effects 
such as increased thermal gradient and bending. 

• Vertical channel placement has only a slight impact on temperature distribution but can have 
a noticeable influence on bending and local induced stresses, as two competing effects are at 
work: i) through thickness temperature gradient and ii) plate asymmetry. A vertical location 
can be determined that will minimize overall bending. 

• Combing both horizontal and vertical channel location movement allows one to balance the 
objective of maximizing cooling while minimizing all detrimental effects. 

Furthermore, two important factors from a modeling viewpoint were also examined in this paper; 
1) the importance of imposing realistic thermal boundary conditions (particularly when evaluating the 
effectiveness of internal cooling) and 2) the degree of required grid refinement. With regard to thermal 
boundary conditions, it was clearly shown, as expected, that if a prescribed surface temperature were 
imposed instead of convective boundary conditions to model the effect of the flame on the plate’s 
boundary, the impact of cooling channels on reducing surface temperatures under the flame was 
artificially reduced. Consequently, an unrealistically high temperature field throughout the plate (which 
also leads to greater predicted stresses and displacements) would result. Similarly it was shown that the 
higher order theory does indeed have some sensitivity toward grid refinement, however this mesh 
sensitivity is not as severe as one might find when using the traditional finite element method. 

Finally, due to the new generalized internal and external boundary capabilities, HOTFGM can 
now be applied to model structures with curved surfaces. This is a significant advance as previously, the 
Cartesian version of HOTFGM was restricted to simulating the response of flat plates. As an illustration 
of this new capability, a cylinder subjected to internal and external pressure and convective thermal 
boundary conditions was modeled and the results compared to an analytical solution. The results 
indicated that HOTFGM can be used to model the thermomechanical behavior of structures with curved 
surfaces, however, some perturbations should be expected in the mechanical results in regions near the 
curved boundaries. 
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Appendix: Thermal/Mechanical Analytical Solution for Thick-Walled Cylinder 

The analytical solution for the cylinder depicted in Fig. 43a is constructed in two steps. First, the 
thermal problem is solved, and then knowing the resulting temperature field distribution throughout the 
cylinder, the mechanical problem is solved. The problem is axisymmetric, and the cylinder is assumed to 
be infinitely long in the out of plane ( z -coordinate) direction. 


Thermal Analysis 

To begin the thermal analysis, consider Fourier’s Law of heat conduction, 


k mr 

d / 


(A. 1 ) 


where q r is the radial component of heat flux per unit area at a point in the cylinder, k is the thermal 
conductivity of the cylinder material, and T(r) is the radially-dependent temperature. No other heat flux 
component is present as the problem is axisymmetric (no variation of temperature in the circumferential 
direction) and the temperature is assumed constant along the length of the cylinder. To obtain the total 
radial heat flow at a particular radial location, Q r , we simply multiply q r by the area over which the heat 
flux occurs; 

Q r = —2nrlk (A. 2) 

dr 

where / is the length of the cylinder. The governing differential equation for the thermal problem is 
obtained by considering the differential volume element shown in Fig. 43b; wherein since the problem 
involves steady-state heat flow, the total heat flow in and out of the differential element must be balanced, 
that is: 


q:=q: ,! (a.3) 

The total heat flow in and out of the differential element can be written in terms of the total heat flow at 
the element midpoint and the radial rate of change of the total heat flow. 


Q‘:=Qr- 


dQ r dr 


dr 2 

Substituting these expressions into eqn. (A.3) yields, 

dQ 

d r 


Qr=Qr + 


= 0 


dQ. r/r 

dr 2 


(A.4) 


(A. 5) 


Substituting eqn. (A. 2) into eqn. (A. 5) yields the governing differential equation for the thermal problem: 


d r 


dT(r) 

d r 


= 0 


(A.6) 


Upon integration we obtain the radial temperature distribution, 

r(r) = c,ln(r) + c 2 (A.7) 

The thermal boundary conditions for the problem are convective at both the inner and outer surfaces 
of the cylinder. The general equation for this type of boundary condition is. 


<7 = h{T-T_) 


(A.8) 
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where h is the convection coefficient of the fluid beyond the boundary, T s is the surface temperature at 
the boundary, and is the temperature of the fluid beyond the boundary. Equation (A.8) is applied to the 
inner and outer boundaries of the cylinder such that the total heat flow at the inner and outer boundaries 
are: 


Q,(r = r l ) = -2xr l h,l[T(r,)-T l ] (A.9) 

Q r (r = r,) = 2xr„hj[T(r„)-T„] (A.10) 


where the subscripts / and o refer to inner and outer, respectively. Note that the sign of the boundary 
condition at the inner surface, (A.9), is reversed as the normal to this surface is negative. Utilizing 
eqn. (A.2) and (A.7) in combination with eqns. (A.9) and (A. 10) the two constants of integration can be 
obtained; they are. 


T-T 


H r i/ r o)~ k {V r i h i + l / r o h o) 


T-T 


c, = 


H^/r^-kil/rfy+l/r^o) 


\ r ; h , 


-Hn) 


+ T 


(A.ll) 
(A. 12) 


Mechanical Analysis 

Assuming elastic, isotropic material behavior the corresponding constitutive equations (in 
cylindrical coordinates) are: 


<*r 

F 

"(1 — v) V V 

£ r -aAT(r) 

Off 

lL 

v (1 — v) V 

£ 0 -ocAT(r) 

o. 

~ (l + v)(l-2v) 

V v (1 — v) 

e° -a AT(r) 


where, £, are the total strain components, <7 are the stress components, A T(r ) is the temperature change 
from a reference temperature, and a . E. and V are the classic elastic material properties (i.e., coefficient 
of thermal expansion, elastic modulus, and Poisson ratio, respectively). The problem is treated as 
generalized plane strain, where £° is the known uniform longitudinal strain throughout the cylinder. 
Consequently, 


a, = Q t e r + Q 2 £ d + Qi £ ° z ~ T AT(r) 
<?e=Q2 £ r + Q\ £ e + Q2 £ °z- r AT ( r ) 

where, 

_E(l-v)_ Ev r = _Ea_ 

^ (1 + v)(l-2v)’ 2 (1 + v)(l -2v)’ (1 -2v) 


(A. 14) 
(A. 15) 


(A. 16) 


The known cylindrical kinematic equations are: 

du u 




(A. 17) 


where u and w are the radial and longitudinal components of displacement, respectively. Since the 
longitudinal strain is known, the last equation in (A. 17) is not needed. Combining eqns. (A. 14) - (A. 17) 
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and substituting the resulting expressions for the radial and tangential stress components into the equation 
of equilibrium. 


d °r ! _Q 

dr r 

yields the following governing differential equation for the mechanical problem: 

d'u 1 du u Yd 

— H ; 

dr r dr r~ Q [ dr 

Equation (A. 19) can then be rearranged to yield. 


[A7X'-)] = 0 


d_ 
d r 


1 d(ru) 


= — — [Ar(r)l 
Q, dr 1 ' J 


r dr 

which can be integrated directly. The resulting expression for the radial displacement is, 

r i 


u(r) = — -- f AT(r) rdr + c^r + — 
Q r { r 


(A. 18) 


(A. 19) 


(A. 20) 


(A. 2 1 ) 


Note that the integral appearing in eqn. ( A . 2 1 ) can be evaluated given the solution to the thermal problem, 
eqn. (A. 7), and the fact that AT(r) = T(r) — T ref , where T rc , is a reference temperature. Evaluating the 
integral, and assigning a function, fj(r), to represent the integral yields. 


r 2 

F t (r) = jAT(r)rdr= 


c i ln(r)-^ + c 2 -T 


ref 


F_ 

2 


c MF)-^ + c 2-T re 


(A.22) 


In order to determine the constants of integration appearing in eqn. (A.21), the mechanical 
boundary conditions must be applied to the cylinder, that is 

oA^r^-P.,, a r (r = r o ) = -P o (A.23) 

Consequently, substituting eqn. (A.21) into eqn. (A. 17) and the results into eqn. (A.14) followed by eqn. 
(A.23) the two constants of integration can be determined and they are: 


C 4 = 


(a-ei-l-4' 

U' rr 

-Q 2 e° z -P i + £±(Q i -Q 2 ) 


(A. 24) 


Q { + Q, 


(A. 25) 


Finally, the uniform longitudinal strain, £° z , can be determined by imposing the constraint that the 
average axial stress, G . , must vanish, that is: 

2nr a 

(A. 26) 


j ° 

f f 0.{r) rdrdO = 0 

Wo- r r)ii 
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The expression for G _(r) is obtained from the third of eqns. (A. 13) wherein E } and are eliminated by 
substituting into eqn. (A. 17) the solution for the radial displacement given in eqn. (A.21). The radial 
integration in eqn (A. 26) can then be performed, and the resulting equation can be solved foi This 
procedure yields. 


£. =■ 


2v 




2/M 

.Mc-'flL 


-2r v-Ea+^-(Q, +Q 2 ) 


(A. 27) 


Thus, given the thermal and mechanical boundary conditions, the uniform £ that results in the vanishing 
of the average longitudinal stress is known. 


N AS A/TM— 2001-2 10702 


28 


Table 1: Material Properties. 


Material 

T(°F) 

E (psi) 

V 

a (10' 6 /°F) 

K (BTU/in s °F) 

Silicon 

77 

4.35x1 0 7 

0.22 

1.83 

4.01x10' 4 

Nitride [10] 

2552 

3.63x1 0 7 

0.19 

1.83 

1.61x1 O' 4 

Mullite [11] 

77 

2.10x10 7 

0.20 

2.94 

7.84x1 O' 5 


2552 

2.10x10 7 

0.20 

2.94 

5.04x10' 5 

Porous 

77 

3.63x1 0 2 

0.25 

6.25 

2.68x1 O' 6 

Zirconia [12] 

2552 

3.63x1 0 2 

0.25 

6.25 

6.69x1 O' 6 


Table 2: Execution times (in seconds) for the different linear equation solvers. 


n b 

Ni 

Number of Subcells 

LEQT1B 

Y12MAF 

uvss 

8 

24 

192 

140 

0.74 

0.44 

16 

48 

768 

3,162 

21.0 

2.07 

26 

74 

1,924 

16,149 

424 

8.51 

48 

136 

6,528 

587,823 



69.23 

96 

272 

26,112 

— 

— 

652 


Table 3: Cylinder parameters and dimensions. 


Parameter 

Value 

r o 

1 in. 

r 

i 

0.8 in. 

T 0 

3600 °F 

T; 

70 °F 

Tref 

70 °F 

K 

0.0003 BTU/in. 2 s °F 

h, 

0.0385 BTU/in. 2 s °F 

P a 

1 000 psi 

P, 

1 0,000 psi 
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20 H ' ■ 1 
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Time to Failure (s) 

Figure 1 : Sample applied stress vs. life curves for Si 3 N 4 [1]. 
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*3 


*1 


Figure 2: 



A geometric model of a heterogeneous graded material in the x 2 - x 3 plane illustrating the 
volume discretization employed in HOTFGM-2D. 
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Generic internal boundary cell 


Figure 3: Illustration of the generic internal and external boundary cells given a functionally graded 
(region L-, by H) internally cooled plate. 
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Fig. 5: Cooled plate dimensions. 
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Figure 6: Cooled plate boundary conditions. Note: units for h are BTU/(in 2 s °F). 
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Figure 7: Five levels of subcell grid refinement: (a) 8x24, (b) 16x48, (c) 26x74, (d) 48x136, (e) 96x272. 
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Figure 8: Temperature fields for the different levels of grid refinement, (a) 8x24, (b) 16x48, (c) 26x74, (d) 
48x136, (e) 96x272. 


N ASA/TM— 200 1 -2 10702 


36 





































0.14 


(a) 



Temperature (°F) 


2800 



0.00 0.10 0.20 0.30 0.40 0.50 

X3 (in.) 

Figure 13: (a) Through-thickness temperature under flame, b) Temperature along top and bottom of 
plate. Effect of grid refinement is shown. 
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Figure 14: (a) Through-thickness x 3 -direction displacement at free edge, (b) x 2 -direction displacement 
along top and bottom of plate. Effect of grid refinement is shown. 
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Figure 15: (a) <X n in substrate along substrate/bondcoat interface, (b) <7,, in substrate along substrate/bondcoat interface, (c) cr„ in substrate 
along substrate/bondcoat interface, (d) in substrate along substrate/bondcoat interface. Effect of grid refinement is shown. 
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Figure 16: Model execution times as a function of the number of subcells in the analyzed geometry 
utilizing three different linear equation solvers. 



(b) 


Figure 17: (a) Subcell grid for the plate with no cooling channels, (b) Subcell grid for the plate with the 
fixed temperature flame boundary conditions. 
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Figure 18: Temperature profiles, (a) No Internal cooling: convective flame boundary condition, (b) Internal 
cooling: convective flame boundary condition, (c) Internal cooling: fixed temperature flame 
boundary condition. 
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Figure 19: I x stress invariant profiles: (a) No Internal cooling: convective flame boundary condition, (b) 
Internal cooling: convective flame boundary condition, (c) Internal cooling: fixed temperature 
flame boundary condition. 
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Figure 20: J 2 stress invariant profiles: (a) No internal cooling: convective flame boundary condition, (b) 

Internal cooling: convective flame boundary condition, (c) Internal cooling: fixed temperature 
flame boundary condition. 
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Figure 21 : Detailed temperature field for the plate with no cooling (reduced range of scale). 
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Figure 22: (a) Through-thickness temperature under flame, (b) Temperature along top and bottom of 
plate. Effects of cooling channel presence and flame boundary condition are illustrated. 


NAS A/TM— 200 1 -2 1 0702 


47 



o 

o 

o 



O 


o 

m 


o 

m- 

o 



cri 

0 

E 

0 

LL 

CO 

0 

O.E 


d 

O) 

C 

c 

c 

0 

X 

CO 

X 

o 

03 

CQ 

0 

E 

o 


O 

CM 

0 

o 

O 

Q. 

E 

03 

c 

d 

c 

o 

LL 

0 

0 



o 

CL 

C 

H 

o 

o 


o 

E 

0 

CO 

~o 

0 

O 

o 

0 

c 

0 

1- 

0 

QQ 

X 

LL 

o 

2 

T — 

o 

"0 

CO 

■O 

0 

1 

i 



0 

X 

T 

T 



CO 

Ll 




o 


I 


_(0 

CD 

C 

c 

03 

x 

O 


o 

o 

O 

o 



o 

in 


o 

M- 


o , 

CO 


CO 

X 


o 

CM 


o 

d 


o 

o 


o 

o ° 

o 

O 

O 

o 

o 

o 

O 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o_ 

o 

o 

o 

o 

o' 

d 

in 

o' 

in 

d 

LO 

d 

in 


in 


M- 

CO 

CO 

CM 

CM 



o 
o 
o 
c> in 


o 

o 

o 


(isd) zzO 


(jsd) «D 



0 

to 

u_ 

00 

-D 

3 

CO 


■O 

c 

CD 

Q) 

O 

c 

CD 

CO 

CD 


CD 

O 

CD 

0 


CD 

o 

o 

■D 

c 

o 

x 


CD 

c 

c 

CD 


O) 

c 

o 

o 

o 


o 

LU 


CO 

XI 

3 

CO 

oo 

c 

o 

CD 

03 

0 

L_ 

00 

X 

=3 

CO 


CD 

O 

0 


CD 

O 

o 

TO 

c 

o 

X 
0 
to 

k_ 

to 

XI 
13 
CO 


b [? 
^ o 

X 


0 

O 

'XI 

0 


0 
0 
to 

4— 

CO 
X 

.£ 5 

15 •£ 

0 r-, 

U «N 

1 0 

I S 
0 
■*-* 

0 
u_ 

to 
X 
3 
CO 

O) 

c 
o 
0 
0 

0 
k_ 

to 
X 
3 
CO 


0 

O 

0 

t: 

0 


0 

o 

o 

“O 

c 

o 

X 

0 
-t— • 
0 
L_ 

to 

X 

3 

CO 


b g> 


CO 

CM 


3 

03 


flame boundary condition are illustrated. 
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Figure 24: (a) Through-thickness x 3 -direction displacement at the free edge, (b) x 2 -direction 
displacement along the top and bottom of the plate. Effects of cooling channel presence and 
flame boundary condition are illustrated. 
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Figure 25: Temperature fields for plates with cooling channels: (a) in the baseline configuration, (b) 
shifted left, (c) shifted right, and (d) with linear spacing. 
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Figure 26: /, stress invariant fields for plates with cooling channels: (a) in the baseline configuration, (b) 
shifted left, (c) shifted right, and (d) with linear spacing. 
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Figure 28. (a) Through-thickness temperature under flame, (b) Temperature along top and bottom of plate, (c) Comparison of temperatures at 
various locations in plate. Effect of shifting the channels horizontally is illustrated. 
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Figure 30: (a) Through-thickness x 3 -direction displacement at free edge of plate, (b) x 2 -direction 
displacement along top and bottom of plate. Effect of shifting the channels horizontally is 
illustrated. 
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(d) ■ ■ ■ ■ ■ ■ 

Figure 31: Plate temperature fields for: (a) baseline configuration, (b) cooling channels shifted up, (c) 
cooling channels shifted down, and (d) minimum bending vertical cooling channel location. 
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Figure 34: (a) Through-thickness temperature under flame, (b) Temperature along top and bottom of plate, (c) Comparison of temperature at 
locations within plate. Effect of shifting the channels vertically is illustrated. 
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Figure 38: (a) Through-thickness temperature under flame, (b) Temperature along top and bottom of 
plate. Effect of horizontal and vertical channel arrangement is illustrated. 
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Fiqure 39' Plate /, stress invariant fields for (a) baseline configuration, (b) cooling channels shifted up 
and to the left, and (c) cooling channels arranged both horizontally and vertically. 


NASA/TM— 2001-210702 


64 












late J 2 s 
id to the I 


i.^i 








C 

o 


=3 c£ 
CD CD 


C 

o 

O 

0 

_C 

0 

CO 

03 

CD 


o& 

Cl 

Z 3 

“O 

0 

> 

O 


♦ * 


o 

in 


o 


u 

CO 

X 

o 

CM 


O 

O 


O 

O 


O 

o 

o 


o 

o 

o 


o 

o 

o 

o 


o 

o 

o 


o 

o 

CD 

cT 

CM 


o 

o 

o 

in 

CM 


o 

o 

o_ 

O 

CO 


o 

o 

o_ 

in 

CO 


o 

o 

o 

o 


(isd) zzo 


c 

0 



o 

LO 


o 

o 


o 

CO ' 


o 

CM 


o 

o 


o 

o 

o 


"A 


c 

0 

E 

0 

o> 

c 

03 


03 

o 
■£ 
0 
> 

- M- «* 

f 3 |> 

O o 5 o 

° o- a 

0) =3 CO 

J= 

0 
CO 
0 
CD 


c 
o 

0 

5. 


u k- 
0 0 
> 0 
o c 




o 

o 

<o 

in 

CM 


o 

o 

o 

o 

CM 


o 

o 

o 


o o 

o o 

o p. 

o LO 


(;sd) 


CO 

X 


o 

m 


o 

d 


0 
4— 1 
0 
k_ 

0 

X 

CO 


o , 

CO 


o 

CM 


o 

o 


o 

o 


o 

o 

o_ 

in 



o 

in 


o 

d 


o 

CO 


0 

o 

t5 

0 

> 

■o 

c 

0 


b S 


c 

o 

N 


0 

o 

0 

Tz 

0 


c 


CO 

X 


o 
0 

§ W 

o 

“O 

c 
o 

-Q 


0 

0 

L. 

■4— > 

CO 

X 

3 

CO 

O) 

c 

o 

0 

0 

*0 

L- 

to 

X 

3 

i o 


0 

o 

0 

•t 

0 


0 

o 

o 

T 3 

c 

o 

X 

0 

0 

l— 

to 

XI 
3 
CO 

CD 

c 

o 

0 


b S 
0 


0 

o 

0 

t: 

0 


CO 

X 

3 

(0 


- b 


CO 

X 


o 

CM 


o 

o 


0 

o 

o 

-O 

c 

o 

X 

0 
4 — • 
0 
k_ 

to 

X 

3 

CO 

O) 

c 

o 

0 

0 

0 

l— 

4-» 

(0 

X 

3 

CO 


b 

"0 


“O 

c 

0 

0 

CO 

t js 
0 0 
f— 4— > 

.!= co 
, 3 


CO 


c 
0 

E 

Q) 0 


0 
k— 
4— 1 

CO 

X 

3 

0 


c £ 
O ^ 
to o 


(isd) ccQ 


-3- 

0 

k_ 

3 

CD 


NASA/TM — 200 1 -210702 


66 



0.14 


(a) 


0.12 

0.10 

-r 0.08 

c 

CM 

X 0.06 

0.04 

0.02 

0.00 

0.0016 



0.0008 

Top 



X3 (in.) 

Figure 42: (a) Through-thickness x 3 -direction displacement at free edge of plate, (b) x 2 -direction 

displacement along top and bottom of plate. Effect of horizontal and vertical channel 
arrangement is illustrated. 
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Figure 45: HOTFGM simulated temperature field in cylinder. 
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Figure 46: Comparison of the HOTFGM simulation and the analytical solution of the through-thickness 
temperature distribution in the cylinder . 
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cylinder simulated by HOTFGM. 
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Figure 48: Comparison of the HOTFGM simulation and the analytical solution for through-thickness 
temperature distribution in the cylinder. 
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Figure 49: HOTFGM simulation of an airfoil cross-section, (a) Subcell grid, (b) simulated temperature field 
with no internal cooling, and (c) simulated temperature field with internal cooling. 
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